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TIMOSHENKO  BEAM  MODELS 


H.T.  Banks  and  J.M.  Crowley 


ABSTRACT 


We  present  cubic  and  linear  spline-based  approximation  schemes  for 
models  of  beams  based  on  the  Timoshenko  theory.  The  schemes  are  used  in 
parameter  estimation  algorithms;  convergence  results  and  numerical  findings 
are  reported. 


PARAMETER  ESTIMATION  IN  TIMOSHENKO  BEAM  MODELS 


Introduction 

When  modeling  large  complex  structures,  it  is  often  desirable  to  model 
them  as  a  simple  structure  such  as  a  beam  or  plate.  For  example,  a  large 
latticed  structure  with  many  components  may  be  efficiently  modeled  as  a 
single  beam  using  either  the  Euler-Bernoulli  or  Timoshenko  theory  (see 
[1] ,  [2],  [3]),  depending  on  the  importance  of  shear  effects  and  rotary 
inertia.  The  identification  of  structural  parameters  (e.g.  bending  and 
shear  rigidity,  rotatory  inertia)  plays  an  important  role  when  modeling 
large  complex  structures  in  this  manner.  Identification  or  estimation  is 
important  in  model  development  and  verification  and  can  also  play  a  funda¬ 
mental  role  in  procedures  for  design  of  feedback  controllers. 

The  purpose  oi  .his  paper  is  to  propose  numerical  methods  for  esti¬ 
mation  of  constant  parameters  in  the  Timoshenko  equations  for  transverse 
vibrations  of  a  beam,  to  report  convergence  results  for  these  methods  and 
to  give  some  numerical  results  demonstrating  the  effectiveness  of  the 
associated  estimation  algorithms.  We  restrict  our  attention  to  the  con¬ 
vergence  properties  of  the  methods  and  do  not  address  such  important  ques¬ 
tions  as  identifiability,  observability,  or  global  stability  of  the  solu¬ 
tions  of  the  underlying  partial  differential  equations.  For  the  presenta¬ 
tion  here,  we  employ  observations  consisting  of  discrete  measurements  of 
transverse  displacements,  although  our  theory  is  sufficiently  general  to 
include  problems  with  measurements  of  velocity  or  strain  in  place  of 
displacements. 


The  theory  and  methods  presented  here  have  also  been  successfully  ap¬ 
plied  to  estimation  of  parameters  in  the  Euler-Bemoulli  equation  [4]  ,  [5] 
and  in  other  hyperbolic  and  parabolic  partial  differential  equations  [6]  - 


[10].  The  equation  to  which  we  restrict  our  attention  here,  however,  is 
the  Timoshenko  system  [11,  p.  300]  : 


(1) 


ytt  =  a(yxx  -  *xJ  +  f(t’x;^ 


'j'tt =  b(yx  ■  w  * 


XX 


0  <  x  <  1,  t  >  0, 


where  y(t,x)  is  the  transverse  displacement  and  tj/(t,x)  is  the  angle  of 
rotation  of  a  cross-section  of  the  beam.  Here  a  =  k'AG/m,  c  =  EA/m, 
b  =  Aa/I  with  A  =  cross  sectional  area,  E  =  Young's  modulus,  G  =  shear 
modulus,  I  *  moment  of  inertia,  k'  =  shear  coefficient,  m  =  mass  per  unit 
lengvU,  and  f  represents  the  transverse  loading  (possibly  depending  on 
unknown  parameters  q) . 

Various  boundary  conditions  are  of  interest,  depending  on  the  physi¬ 
cal  situation.  Among  these  are: 


(a)  y(t,0)  *  <Kt,0)  =  y(t,l)  =  *Kt,l)  =  0  (fixed  end  beam) 

(2)  (b)  y(t,0)  =  ij)  (t,0)  =  y(t,l)  =  ip  (t,l)  =  0  (simply  supported  beam) 

(c)  y(t,l)  a  Kt.l)  =  y„(t,0)  -  <Kt,0)  =  (t,0)  =  0  (centilevered 

beam) . 


In  addition  one  must  specify  initial  conditions  y(0,x),  yt(0,x),  *K0,x), 
*t(0,x)  for  (1). 

Suppose  that  we  are  given  a  set  of  observations  (y_)  that  corres¬ 
pond  to  measurements  of  the  transverse  displacements  (y(t^,Xj)}  at  times 
t^  and  at  points  x^ .  We  assume  that  the  structure  that,  we  are  observing 
can  be  modeled  by  equation  (1).  We  seek  those  values  of  certain  unknown 


3 


parameters  q  (either  unknown  coefficients  in  (1)  and  unknown  parameters 
in  the  load  function  and/or  in  the  initial  conditions)  so  that  the  solution 
(y,<p)  =  (y(q)  .’/'(q) )  of  (1)  best  fits  the  data  (y^)  in  some  sense.  Here 
we  shall  measure  the  fit  of  the  model  to  the  data  by  a  sum  of  squares 
functional  of  the  form: 

J(q.y)  =  l  IPji  -  y(vV I2, 

i  J 

The  parameter  identification  problem  is  then  to  find  from  a  given  parameter 
set  Q,  values  q  of  the  parameters  which  minimize  J(q,y)  subject  to  y 
being  a  solution  of  (1) . 

Since  explicit  solutions  of  (1)  are,  in  general,  impossible  to  find, 
we  are  led  to  seek  approximate  solutions  of  the  parameter  estimation  prob¬ 
lem.  We  first  approximate  solutions  of  (1)  by  some  appropriate  method  to 

N  N 

obtain  approximate  solutions  (y  ,\p  ).  The  associated  estimation  problem 

—N 

then  consists  of  finding  a  vector  of  parameters  q  which  minimizes 
N 

J(q,y  )•  Our  goal  here  is  to  present  a  class  of  methods  for  estimation  of 

parameters  for  which  we  have  proven  convergence  of  the  solutions  of  the 

— N  — 

approximating  problems.  That  is,  we  have  shown  that  q  — ^-q  as  N  ». 
Further,  we  hope  to  demonstrate  numerical  feasibility  of  these  methods. 

The  class  of  methods  we  examine  for  approximating  solutions  to  (1) 
are  projection  or  Ritz-Galerkin  type  methods.  The  initial -boundary  value 
problem  for  (1)  is  rewritten  in  terms  of  a  system  of  first  order  operator 
differential  equations  in  an  appropriate  Hilbert  space  Z.  This  abstract 
system 

z(t)  =  A(q)z(t)  +  F(q,t) 

(3) 

2(°)  =  zQ 


4 


is  then  approximated  by  a  system  of  the  form 


'(t)  =  AN’(q)zN(t)  ♦  Aq.t) 


z  (0)  =  P  zQ, 


N  N  N  N  N  N 

with  A  (q)  =  P  A(q)P  and  F  =  P  F,  where  P  is  the  orthogonal  pro- 

N 

jector  of  the  Hilbert  space  Z  onto  some  finite  dimensional  subspace  Z  . 

N 

We  shall,  in  each  case  discussed  below,  choose  Z  as  the  linear  span  of 
a  particular  set  of  spline  functions  satisfying  the  appropriate  boundary 
conditions.  For  these  schemes,  we  have  established  convergence  as  N  -*■  <® 
of  the  parameter  estimates  to  their  best-fit  values  (for  details,  see 
[4],  (SI,  (61). 


Formulation  of  the  Abstract  System 

The  Timoshenko  equations  (1)  can  be  rewritten  as  a  first  order  system 
of  partial  differential  equations  in  a  number  of  ways.  Certain  reformu¬ 
lations  are  more  suitable  to  handling  given  boundary  conditions  than 
others.  Each  reformulation  leads  to  a  different  abstract  system  (that  is, 
the  operator  A  in  (3)  and  perhaps  the  Hilbert  space  Z  depend  on  the 
particular  reformulation)  and  thus  can  lead  to  a  different  approximating 
system.  We  examine  two  such  reformulations  here.  The  first  we  investi¬ 
gate  is  most  natural  when  considering  beams  with  fixed  ends  (i.e.  boundary 


conditions  (2a)).  We  formally  rewrite  (1)  with  z  =  (vi»v2>V3>v4)  = 
Cy»yt»'Mt)  as 


5 


JTV1  =  v2 


(5) 


wv2  =  a  7TV]  -  a-kv 3  +  f(t*x;^ 

3x 


at  v3  =  v4 


a  ,  a  a- 

—  V  ,  =  b  -5 -  V ,  +  C  - zr  V.  -  bv.  . 

at  4  a.\  i  Sx2  ^  3 


Letting  z^t)  =  v^t,-).  this  formulation  leads  to  the  abstract  equation 
(3)  in  2  =  hJ  x  L2  *  hJ  x  L2,  with  q  =  (a,b,c,q)  =  (q^q^q^q) ,  where 


(6) 


A(q)  = 


qiD 


l  <2° 


0  I 

2 


-*1* 


0 

0 

I 


q3^-q2I  0 


D  = 


on  Dom(A(q))  =  (H2  n  hJ)  x  l2  x  (H2  fl  hJ)  x  L2,  and 


F(q.t)  = 


f(t,-;q) 

0 

0 


Cantilever  boundary  conditions  cannot  be  handled  easily  in  this  for¬ 
mulation;  another  reformulation  is  more  appropriate.  Using  the  change  of 
variables  Vj  =  yx  -  p,  v2  =  y^,  =  ^x>  v4  =  ^t»  v5  =  Y»  one  can  obtain 
a  system  of  first  order  partial  differential  equations 


6 


9Vj  3v2 

w =  '  v4 

3v_  3v 

~  =  a  33T  +  f(t'x;^ 


(7) 


at 

3V2 

at 

3v4 

at 

at 


3x 

sv 

bvi  *  c  ST 


=  V 


2’ 


which,  under  the  correspondence  z^t)  =  v^t,-).  leads  to  the  abstract 

2  2  2  2  2 

equation  (3)  in  Z=L  *  L  x  L  x  l  x  l  ,  with  q  =  (a,b,c,q)  where 


and 


(8) 


F(q.t)  = 


A(q)  = 


V 


D 

0 

0 

0 

I 


f(t, • ;q) 
o 
o 
o 

o 
o 
0 


q3° 


-I 

0 

D 

0 

0 


0 

0 

0 

0 

0 


on  Don(A(q) )  *  HJ,  x  hr  x  Hl  x  Hr  x  Hr-  Here  the  sPaces  H*  =  {<f>  €  H1: 
♦CO)  *  0}  and  H*  =  {$  €  H*:  $(1)  =  0}  are  defined  to  deal  with  the 
appropriate  boundary  conditions  v^t.O)  =  v3(t,0)  =  v2(t,l)  =  v4(t,l)  = 
vs(t.l)  =  0. 

One  nay  view  the  techniques  employed  to  approximate  the  partial  dif¬ 
ferential  equations  above  as  a  particular  Galerkin  approximation  to  (5)  or 


t 


(7).  Our  formulation  as  an  abstract  equation  in  a  Hilbert  space  Z 
provides  a  convenient  setting  in  which  to  prove  convergence  of  the  systems. 


Convergence  Results 

As  we  have  already  stated,  once  we  have  the  initial -boundary  value 

problem  for  (1)  rewritten  as  an  abstract  system  (3)  in  an  appropriate 

N  N 

Hilbert  space  Z,  we  approximate  (3)  using  projections  P‘  :  Z  — >Z'  .  In 

N 

particular  we  choose  a  sequence  of  finite  dimensional  subspaces  Z  ,  with 
N  N 

Z  c  Dom(A),  such  that  P‘  z  ■+•  z  for  all  z  6  Z,  and  use  the  sequence  of 
approximating  ordinary  differential  systems  (4). 

Given  this  approximation,  we  can  formulate  an  approximate  identifica- 

♦  •  “"W  ^ 

tion  or  estimation  problem:  find  q‘  which  minimizes  J(q,y  )  over  Q, 

N  N  M 

where  y  (t,-)  =  Zj(t)  with  z  the  solution  of  (4). 

We  outline  briefly  (for  further  details  on  the  theory,  see  [4],  [5], 

[6],  [12])  the  essential  ideas  involved  in  establishing  convergence  as 
N  -*■  ®  of  the  estimates  q^’  to  best-f’t  parameter  values  q  for  (3). 

In  our  approach,  the  concept  of  dissipativeness  plays  a  central  role.  An 
operator  A  is  said  to  be  dissipative  ([13],  [14],  [15])  if  there  exists  a 
constant  w  such  that  <Az,z>  <_u)<z,z>  for  every  z  in  Dom(A) .  It  is 
this  property  of  the  operator  A(q)  (and  not  self-adjointness)  which  we 
require  when  formulating  the  abstract  system  and  making  convergence  argu¬ 
ments. 

To  prove  convergence  one  first  shows  that  A(q)  satisfies  a  uniform 
(in  q)  dissipative  inequality  in  Z  and  that  A(q)  generates  a  CQ 
semigroup  T(t;q).  In  the  cases  of  the  operators  in  (6)  and  (8),  the  latter 
is  accomplished  using  standard  results  on  perturbation  of  semigroups.  The 
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N  N  N  N 

approximating  operators  A  (q) ,  which  are  defined  by  A  (q)  =  P  A(q)P  , 

are  easily  shown  to  be  uniformly  (in  N  and  q)  dissipative  and  to  gen- 

N 

At 

erate  a  stable  family  of  approximating  solution  semigroups  e  such 
that  |  expfA^CqJt]  j  <.Me<1>t,  where  M  and  w  are  independent  of  N  and 
q  for  q  in  a  bounded  set.  Standard  estimates  from  spline  interpolation 
theory  may  be  used  to  obtain  bounds  on  the  error  of  the  projections 
|PN'z  -  z|.  The  projection  error  bounds  then  in  turn  are  used  to  show 
that  AN(q)z  — »  A(q)z  in  an  appropriate  sense.  One  then  employs  the 
Trotter-Kato  theorem  [13,  p.  90]  (a  functional  analytic  version  of  the  Lax 

Equivalence  Theorem;  stability  +  consistency  convergence)  to  argue  that 

N  N 

[exp  A  (q)t]  z  — >T(t;q)z  and,  moreover,  that  z  (t;q)  — z(t;q)  where 

N 

z  and  z  are  solutions  of  (4)  and  (3)  respectively.  These  convergence 

— N 

results  can  then  be  used  to  argue  that  any  sequence  q  of  solutions  to 
the  approximating  estimation  problems  (involving  (4))  has  a  subsequence 
q^k  that  converges  to  some  q  that  is  a  solution  of  the  original  esti¬ 
mation  problem  for  (3) . 

Implementation  and  Numerical  Results 

As  noted  above,  to  implement  our  methods  one  begins  by  writing  the 
underlying  initial -boundary  value  problem  as  an  abstract  system  in  an  ap¬ 
propriately  chosen  (Hilbert)  state  space.  To  do  this,  one  writes  the 
partial  differential  equation  (1)  in  the  form  of  a  system  of  first  order 
partial  differential  equations  (such  as  (5)  or  (7))  and  chooses  the  ap¬ 
propriate  state  space  in  which  to  pose  the  problem  as  an  abstract  initial 
value  problem.  Of  course,  one  must  check  that  the  chosen  formulation  leads 
to  a  well-posed  problem  and  that  the  resulting  operator  A(q)  generates  a 


*7 
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CQ  semigroup  T(t;q)  in 
in  this  paper). 


(we  have  done  this  for  all  cases  discussed 


Given  this  formulation  in  a  space  Z,  one  then  chooses  the  approximat- 
N 

ing  subspaces  Z  in  Dom(A(q)).  A  standard  Galerkin  approach  may  be 

used  to  obtain  a  concrete  realization  of  the  approximating  system  (4). 

12  12 

To  illustrate,  we  first  consider  the  equation  (5)  in  Z  =  *  L  *  *  L  . 

We  use  as  basis  for  our  approximating  space  Z  a  set  of  cubic  splines 
which  satisfy  the  boundary  conditions. 


To  be  more  precise,  given  a  partition  A 


N  {xi}?=o’  xi  =  i/N’  of 


[0,1],  a  cubic  spline  is  any  function  s  €  [0,1]  such  that  s  is  a 

cubic  polynomial  on  each  subinterval  (x^,xi+1),  i  =  0,...,N-1.  We  denote 
by  S^(AN)  the  set  of  cubic  splines  with  knots  AN.  A  cubic  spline  inter- 
polant  to  a  function  f  €  C  on  [0,1]  is  the  function  s  €  S  (A  )  satis¬ 
fying  Ds(n)  *  Df(n)  for  p  =  0  and  1,  and  s(Xj)  =  f(Xj),  i  =  0,1,..., N. 
3  N 

As  a  basis  for  S  (A  ) ,  the  standard  B-splines  [16,  p.  89]  can  be 
employed:  Define 


C (x)  =  l  (-1)  02)(*-i >  +  .  x  €  [-2,2] 

i=-2  1 


where 


Then  define 

(9)  C?(x)  =  C((x-x.)/h) 

where  h  *  1/N.  Defining  an  extended  partition  A*  =  the  set  of 

cubic  B-splines  provide  a  basis  for  S^(AN).  In  anticipation  of 

our  consideration  of  boundary  conditions  of  type  (2a)  in  our  ensuing  dis- 


10 


cussions,  we  define  S3(AN)  =  {s  £  S3(A^)|s(0)  =  s(l)  =  0}.  As  a  basis  for 

S3(AN)  we  take  {CN}N  n  where 
0  1  x=0 


N  =  0N 

Li  V 


„N 


4C 


2  <_  i  <_  N-2, 

;N 

-1 


c„  =  c„  -  4cf!  . 

N  N  N+l 


_N  -N 
C1  =  C0 


■'N 

4C1 


K  *N 

C  =  C  -  4C 
N-l  N  N-l 


For  the  approximating  subspaces  Z^  we  choose  ZN  =  S3(AN)  *  S3(A^)  * 
Sq(AN)  x  Sq(AN).  With  C1^  as  the  basis  elements  just  defined  for  Sq(AN), 
we  obtain  as  a  basis  for  ZN  the  set  {6^}4^q3  where 


(C^, 0,0,0),  i  =  0.....N 

N  J  (0'Ci-(N*ir°’0)’  1  =  N+l, . . .  ,2N+1 

1  j  (0.0,C^_(2N+2),0),  i  =  2N+2, . , . ,  3N+2 

(0,0,0,C^_(3N+3)),  i  =  3N+3, . . . ,4N+3. 

N  N  N 

Given  Z  ,  in  the  Galerkin  approach  to  our  method  one  seeks  z  €  Z 

satisfying 


(10) 


N  N  N  N  N 

<z  -  A  z  -  P  F,6>  =0  for  all  8  €  Z  , 


where  <,>  is  the  inner  product  in  Z.  it  suffices  to  require  (10)  for 

B  =  j  =  0,  . . . ,  4N+3 .  Since  zN  €  ZN  it  can  be  represented  by 

zN  =  £  w^(t)6^  and  the  condition  (10)  thus  reduces  to 


(ID 


.NaN  r  Na0N  _ 
<)  w.B.  -  L  w.A0.  -  F, 

L  i  l  L  l  l  ’ 


for  j  =  0, 1, . . . ,4N+3.  This  Ritz-Galerkin  argument  leads  to  a  4N+4- 
dimensional  matrix  system  for  the  coefficients  w^(t)  in  the  expansion 

for  zN(t)  relative  to  the  basis  for  ZN.  In  particular,  one  finds  that 

N  N  N 

w  (t)  =  (w0(t),...,w4N+3(t))  satisfies 


with  (Q*1)^  =  <8i>8>,  (1^)  =  <8i,A6j>,  (RNF).  =  <B.,F>.  This  becomes 

(12)  wN(t)  =  GNwN(t)  +  F^ 

with 


where  =  (f^,...,f^)  corresponds  to  the  "load"  f  in  (5)  and  is 

N  N  - 1  /SN  N 

given  by  r  =  (Aj)  (R  f)  with  (R  f)^  =  <CY,f>.  Finally,  since 

z^(t)  -  y(t,-).  the  approximation  to  the  displacement  y  is  given  by 

N 

(13)  yN(t, •)  =  z"(t)  =  l  w. (t)C^ . 

1  i=0  1  1 

We  next  summarize  results  of  some  of  our  numerical  experiments  using 
the  approximation  scheme  based  upon  cubic  splines  just  described.  But 
first  we  describe  how  our  numerical  tests  were  performed.  We  took  as  "data" 
the  values  of  the  solution  of  a  model  of  the  form  (1)  whose  parameters  lq 
were  known  (i.e.  chosen  by  us)  and  sought  a  solution  of  the  approximate 

identification  problem  for  different  values  of  N.  Specifically,  the  "data" 
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{y_},  i  =  l,...,r,  j  =  1,...,£  were  generated  in  this  case  using  the  LZ 

formulation  of  the  Timoshenko  equations  (5)  and  using  a  general  purpose 

computer  code  (M0L1D)  to  solve  this  system  of  first-order  hyperbolic 

equations  to  obtain  displacements  y(t,x).  Then  {y^},  i  =  l,...,r, 

j  =  1,...,£  with  y  =  y(t^,Xj)  were  used  as  the  observations  or  input 

to  the  approximate  identification  package.  For  the  numerical  experiments 

detailed  below,  we  used  £  =  9,  r  =  10  and  generated  the  "observations" 

at  Xj  =  j/10,  j  =  1,...,9  for  times  t^  =  i/10,  i  =  1,...,10. 

For  each  fixed  N,  an  optimization  algorithm  (the  Levenberg-Marguardt 

in  the  IMSL  package  ZXSSQ)  was  employed  to  solve  the  approximate  estima- 

N  0 

tion  problem.  This  optimization  algorithm  requires  an  initial  guess  q  ’ 
for  the  parameters,  which  is  referred  to  as  the  "start  up"  value  in  the 
tables  below.  The  values  of  the  parameters  to  which  the  optimization 

— M  — M 

algorithm  converged  for  a  given  N  are  denoted  by  q  and  J(q  )  is 
the  cost  functional  (residual  sum  of  squares)  for  those  values  of  the  para¬ 


meters,  while  the  values  q  are  called  "true  values". 


We  consider  the  motion  of  a  beam  initially  at  rest  with 


fixed  ends.  That  is,  our  beam  is  described  by  the  system  (1)  with 

-2t 

a  =  q^,  b  =  q£,  c  =  q3,  f  =  lOe  sin  2t,  boundary  condition  (2a)  and 
vanishing  initial  data.  The  numerical  results  are  given  in  Table  1. 


Table  1 


N 

$ 

qN 

3 

6 

.9882 

726.19 

3.6864 

.526 

X 

XJ- 

1 

o 

r-H 

8 

.9969 

781.00 

3.9218 

.735 

X 

10-5 

10 

1.0009 

794.29 

3.9684 

.108 

X 

10-5 

12 

1.00036 

794.90 

3.9732 

.165 

X 

NO 

1 

o 

r-H 

16 

1.00033 

797.85 

3.9883 

.256 

X 

*-» 

o 

1 

true  value 

1.0 

800. 

4.0 

START  UP 

.9 

1000. 

3.9 

Example  2.  We  consider  a  clamped  beam  deformed  to  the  shape 
<|>(x)  =  cos  Ax  +  cosh  Ax  -  K(sin  Ax  +  sinh  Ax)  with  A  *  4.730, 

K  *  (sin  A  +  sinh  A)/(cos  A  ♦  cosh  A),  then  allowed  to  vibrate  freely, 
which  can  be  described  by  the  system  consisting  of  equation  (1)  with 
f  =  0,  initial  conditions  y(0,x)  =  $(x),  yt(0,x)  =  0,  <|>(0,x)  =  4>'(x), 
^(O.x)  =  0,  and  boundary  conditions  (2a).  The  numerical  results  are 
given  in  Table  2.  As  in  Example  1,  the  methods  perform  exceedingly  well 


Table  2 

-N 

-N 

5 

4 

1.1812 

1325.3 

3.6437 

.613  x  i(f2 

8 

.9480 

1125.1 

3.993 

.167  x  10‘2 

16 

.9908 

1152.3 

3.8875 

.242  x  10'3 

24 

1.0009 

1193.7 

3.9771 

.122  x  IQ-3 

TRUE  VALUE 

1.0 

1200. 

4.0 

START  UP 

.9 

800. 

3.8 

14 


We  next  consider  examples  for  the  formulation  (7),  which,  as  we  have 
already  indicated,  is  very  convenient  when  treating  cantilevered  beams. 
Recall  (see  (8))  in  this  case  Dom(A(q))  =  X  X  HR  X 

that  the  cantilever  boundary  conditions  are  Vj(t,0)  =  Vj(t,0)  *  V2(t,l)  = 
v4(t,l)  =  v^(t,l)  =  0.  For  this  formulation  one  can  easily  use  either 
cubic  or  linear  splines.  We  first  outline  the  procedure  using  cubic 
elements. 

3  N 

The  set  of  cubic  splines  S  (A  )  and  its  basis  of  N+3  B-splines 

r^NiN+l 

{CL  ^  were  defined  in  (9)  above.  These  were  then  modified  to  yield 

basis  elements  satisfying  particular  boundary  conditions.  We  do  the 

N 

same  again  here.  To  obtain  a  basis  for  Z  whose  elements  satisfy  the 
cantilever  boundary  conditions  (observe  that  the  projection  method  requires 
that  they  be  in  Dom(A)),  we  define 


S3(AN)  =  (s  €  S3(AN):  s(0)  =  0} 

S3(AN)  =  {s  €  S3(AN):  s (1)  =  0} 

and  as  a  basis  for  sf(AN)  take  (cWi  where  C ^  =  C^,  2  <  i  <  N+l, 

L  l  1=0  i  i  —  —  9 


(the 


CL  are  the  B-splines  given  in  (9)  above). 
As  a  basis  for  S3(AN)  take  {B?}^+*  where 

K  1  1=U 
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N  3  3  3  3  3 

Now  we  take  Z  =  S,  x  S_  x  S  x  sD  x  S_  in  the  case  of  cubic  spline 

L  K  L  K  K 

approximations . 

Turning  to  approximation  schemes  employing  linear  splines,  we  use 
usual  "hat"  functions 

'  N(x  -  xi _1),  x4_j  <  x  <  xt 

4?(x)  »<  N(x.+1  -  x),  x.<x<x.+1 

0  ,  otherwise 

where  AN  =  {x^}^_0  are  the  knots  as  defined  in  setting  up  the  cubic 

splines.  Denote  by  S1^)  the  set  of  linear  splines  with  knots  AN; 

1  N 

i.e.,  s  6  S  (A  )  is  a  function  s  €  C[0,1]  such  that  s  is  a  linear 
function  on  each  subinterval  [x^x^],  i  =  0,...,N-1.  The  N+l  func¬ 
tions  form  a  basis  for  S1(AN).  Define  SJ(AN)  = 

span{£j, . , . and  S*(AN)  =  spaniZ^, . . . j}.  Then  S*  and  S* 
are  N-dimensional  subspaces  of  linear  splines  in  H*  and  H*  respectively. 
We  take  =  S*  x  S*  x  x  s*  x  for  the  linear  spline  approximations. 

The  usual  Ritz-Galerkin  formulation  leads  to  a  5N+10-dimensional 
matrix  system  of  ordinary  differential  equations  for  the  coefficients 
w^(t)  in  the  expansion  for  z^(t)  relative  to  the  cubic  spline  basis 
for  ZN  and  a  5 N-dimensional  system  in  the  linear  spline  case.  For  the 
cubic  spline  case,  the  approximating  matrix  equations  (a  system  of  5N  ♦  10 
ordinary  differential  equations  in  this  case)  again  have  the  form 

QNwN(t)  =  KNwNCt)  ♦  FN{q,t) 

where  now  QN  and  KN  are  5N  ♦  10  square  matrices  defined  by 
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with 


and 


with 


qn  = 


(Q*)..  =  <C?,CN> 

1  ij  1  j 


(Q23ij  =  <Bi*BjS  i>i  *  0 . N+1» 


-V<>' 

0 


K4 


q2(^)T 


0 

0 

0 


0 

0 

0 


.ytt 

*5 


k" 

*4 


„N 

Q2 


0 

0 


(K^l)  .  .  =  <C^,B^> 
3  ij  l  j 


(Kj)ij  *  <C^,DB^>,  i,j  -  0, . . . ,N+1, 


(again,  here  <  ,  >  is  the  inner  product  in  l^) ,  and 


F*  *  (O/.O.O.O)7, 
(f\  =  <f  (t ,  •  ;q) ,  B^>, 


,N 


and  f(t,x;q)  is  the  "load"  term  in  (7).  The  matrices  Qi  are  7- 


banded  and  symmetric  while  the  matrices  K?  are  7-banded  but  not  sym^ 


metric. 
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For  the  linear  spline  case,  the  approximating  system  of  5N  ordinary 

N 

differential  equations  has  the  same  structure  as  above  with  C  replaced 
by  i  =  0.....N-1,  and  by  ij,  i  =  0,...,N-1.  The  5N  x  5N 

matrices  and  however  are  only  tridiagonal  in  the  linear  case. 

Example  3.  In  our  first  example  with  cantilever  boundary  conditions 
presented  here,  we  used  values  for  the  coefficients  reported  in  [1]  where 
a  cantilever  truss  system  is  "modelled"  as  a  continuous  beam.  That  is, 
we  choose  as  true  values  a  =  =  298.257,  b  =  =  413520.,  c  =  = 

3590.57  in  the  equations  (1)  with  f  =  exp(-2t-2x),  initial  conditions 
for  a  beam  at  rest,  and  with  cantilever  boundary  conditions  (2c).  Again, 
simulated  data  (with  observations  at  x..  =  j/10,  t^  =  i/10,  j  =  0,...,9, 
i  =  1, . . . ,10)  were  generated  by  solving  (7)  using  the  general  purpose  PDE 
package  M0L1D. 

Of  course,  in  practice  one  might  be  interested  in  using  fewer  sensors 
or  points  x..  where  data  is  sampled.  Recall  that  our  original  intent  is 
to  examine  the  convergence  properties  of  a  numerical  parameter  estimation 
scheme  and  we  wish  to  avoid  questions  of  observability,  identifiability, 
and  optimal  location  of  sensors.  As  a  comparison,  however,  we  have  in¬ 
cluded  below  another  example  where  data  is  sampled  at  one  point  only 
(the  tip  of  the  cantilever).  The  following  three  tables  summarize  the 
results  for  this  example  using  both  the  cubic  and  linear  spline  schemes. 

In  each  of  these  experiments,  only  one  parameter  was  treated  as  unknown 
and  the  other  two  were  fixed  at  their  true  value. 


Table  3.1 


(treating  q^  as  the  unknown  parameter) 


N 

Linear  Spline 
-N 
•*1 

Cubic  Spline 

-N 

qi 

2 

339.913 

297.792 

3 

283.762 

298.032 

4 

310.843 

298.109 

5 

292.163 

6 

300.690 

8 

302.266 

10 

299.072 

TRUE  VALUE 

298.257 

298.257 

START  UP 

200. 

200. 

Table  3.2 

(treating  q2  as  the  unknown  parameter) 

Linear  Spline  Cubic  Spline 


N 

rN 

q2 

^2 

1 

415363 

2 

351677. 

414597 

3 

413632 

4 

391564. 

413584 

5 

398172. 

6 

405907. 

8 

409217. 

10 

410404. 

TRUE  VALUE 

413520. 

413520, 

START  UP 


360000 


360000 
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Table  3.3 

(treating  as  the  unknown  parameter) 


N 

Linear  Spline 

43 

Cubic  Spline 
qN 

43 

1 

3553.57 

2 

4261.79 

3590.95 

3 

3410.55 

3587.77 

4 

3779.85 

3590.50 

6 

3684.83 

10 

3541.13 

TRUE  VALUE 

3590.57 

3590.57 

START  UP 

3000.00 

3000.00 

We  note  that  both  the  linear  and  cubic  spline  approximations  converge 
rapidly  and  give  reliable  estimates  even  for  small  values  of  N. 

The  optimization  portion  of  the  computer  codes  (IMSL's  ZXSSQ) 
encountered  difficulties  in  this  example  when  we  attempted  to  estimate 
all  three  parameters  simultaneously,  due  to  the  large  difference  in  mag¬ 
nitude  of  the  parameters.  This  can  be  corrected  by  a  rescaling  in  the 
optimization  routine. 

Example  4.  We  consider  next  another  example  involving  cantilever 
boundary  conditions.  For  this  example  we  compare  not  only  cubic  versus 
linear  spline  approximations,  but  also  compare  results  with  a  large  num- 

9 

ber  of  sensors  (ten  points,  x^  =  i/10)  to  those  obtained  in  the 

case  of  a  small  number  of  sensors  (one  point,  x^  =  0). 

The  Timoshenko  equations  (1)  with  the  beam  initially  at  rest  and  with 
cantilever  boundary  conditions  (2c)  was  the  model.  The  applied  load 
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f  had  the  form  f(t,x)  =  exp(-5t)exp(-4x)sin(2t) .  Simulated  data 
was  again  generated  employing  the  general  purpose  PDE  solver  MOL1D  using 
the  "true  values"  a  =  qj  =  1 . 0,  b  =  =  1 . 2  *  103,  c  =  q3  =  12.0.  We 

note  that  parameter  values  on  this  scale  are  obtained  when  a  second  change 
of  variables  (in  addition  to  rescaling  the  spatial  variable  to  the  inter¬ 
val  [0,1]  as  in  [1])  is  performed  by  rescaling  the  time  variable  t  =  at' 

2 

where  t*  is  the  original  time  scale,  and  a  is  of  the  same  order  of 


magnitude  as  q^.  Two  sets  of  "observation"  data  were  employed.  Both 
used  the  value  of  the  approximate  solution  y_  =  y(t^,x^)  at  t.  =  i/10, 
i  =  1,...,10;  the  number  of  spatial  points  used  were: 


(A)  ten  points  x^  =  i/10,  i  =  0,...,9. 

(B)  one  point  xQ  =  0  (at  the  free  tip) . 

The  following  tables  summarize  the  results  obtained.  Table  4.1  through  4.4 
involve  use  of  data  of  the  form  (A),  while  Tables  4.5  through  4.8  result 
from  use  of  data  of  type  (B) .  Thus  Tables  4.5  -  4.8  present  results  using 
only  observations  at  the  tip  (free  end)  of  the  cantilever  beam  at  ten 
times  (tA  =  i/10,  i  =  1,...,10). 

Table  4.1 

(treating  q^  as  the  unknown  parameter) 


Linear  Spline  Cubic  Spline 


N 

-N 

<1 

< 

1 

1.05917 

2 

1.01066 

3 

1.52290 

.99549 

4 

1.18569 

1.00094 

5 

1.15071 

6 

1.05861 

8 

1.02103 

12 

1.00702 

TRUE  VALUE 

1.0 

1.0 

START  UP 

.8 

.8 
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Table  4.2 

(treating  q2  as  the  unknown  parameter) 


Linear  Spline  Cubic  Spline 


1 

1.10840 

2 

1.19366 

3 

.76109 

1.208S3 

4 

.98753 

1.19892 

5 

1.02709 

6 

1.12289 

1.20043 

8 

1.17090 

12 

1.18983 

16 

1.19399 

TRUE  VALUE 

1.2  x  103 

1.2  x  io3 

START  UP 

.8  x  103 

.8  x  IQ3 

Table  4.3 

(treating  q3 

as  the  unknown  parameter) 

Linear  Spline 

Cubic  Spline 

N 

H3 

55 

1 

12.9822 

2 

12.0636 

3 

19.1429 

11.9116 

4 

14.5998 

12.0110 

5 

14.0550 

12.0118 

6 

12.8214 

11.9956 

8 

12.2942 

12 

12.0991 

16 

12.0581 

TRUE  VALUE 

12.0 

12.0 

START  UP 


10.0 


10.0 
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Table  4,4 

(treating  all  three  parameters  as  unknown) 


Linear  Spline 

N 

-N 

<1 

C*  103) 

-N 

^3 

J(F) 

3 

1.68674 

1.61611 

14.4218 

. 197*10  ? 

4 

3.38199 

1.68874 

5.4154 

.166x10", 

6 

1.06937 

.65044 

6.3526 

.478x10", 

8 

.99744 

.78523 

7.9986 

.152x10'^ 

12 

1.00496 

.98184 

9.2938 

.328x10": 

16 

1.00972 

1.03337 

10.2396 

.9r^xl0"b 

TRUE  VALUE 

1.0 

1.2  x  103 

12.0 

START  UP 

.8 

.8  x  103 

10.0 

Cubic  Spline 

N 

-N 

*1 

103) 

$ 

J  (?) 

2 

2.29673 

1.65908 

6.6929 

.561x10'* 

3 

1.27137 

1.16502 

8.6100 

.161x10": 

4 

1.01927 

1.19287 

11.6788 

.507x10", 

6 

1.01317 

1.20155 

11.9933 

.957x10"' 

.470x10" 

8 

1.00067 

1.20160 

12.0051 

TRUE  VALUE 

1.0 

1.2  x  103 

:.2.o 

START  UP 

.8 

.  8  x  103 

10.0 

Table  4.5 

(q1  treated  as  unknown) 

Linear  Spline  Cubic  Spline 


N 

Sf 

J<?) 

1 

1.08757 

.563*10"3 

2 

7 

.97849 

.  194x10"'? 
.377x10": 

3 

1.43945 

.950x10", 

.97627 

4 

1.27313 

.151x10", 

1.00126 

.967xl0_i> 

5 

1.16280 

.744x10", 

6 

1.09265 

. 647x10" l 

8 

1.04641 

.2S6X10"5 

TRUE  VALUE  1.0 


1.0 


START  UP 


.8 


.8 
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Table  4.6 


(q2  treated  as  unknown) 


Linear  Spline 


N 

q^(*103) 

JCq”) 

q^cxio3; 

1 

1.07420 

2 

-3 

1.23856 

3 

.810693 

.108x10  ° 

1.23247 

4 

5 

.919741 

1.01736 

.158X10  , 
.747x10", 

1.19859 

1.19745 

6 

8 

1.08323 

1.14092 

.636x10, 

.255x10", 

1.20103 

10 

1.16499 

.111x10"^ 

12 

1.17608 

.SSlxlO"4 

16 

1.18724 

.180x10 

Cubic  Spline 


T  /“ N\ 

J(q  ) 


.502x10 

.189x10* 

.370x10* 

.967x10* 

.105x10* 

.131x10* 


TRUE  VALUE  1.2x10' 


STARTUP  1.0x10' 


1.2x10' 


1.0x10' 


Table  4.7 


(q3  treated  as  unknown) 


Linear  Spline 


Cubic  Spline 


N 

J(q  ) 

-N 

<3 

.,-Nx 

J(q  ) 

1 

13.3809 

,512*10"3 

2 

3 

17.8588 

.118x10'^ 

11.5821 

11.6669 

.  365x10* l 

4 

5 

15.6401 

14.1666 

.165x10, 

.774x10", 

12.0138 

.968x10" 

£ 

6 

8 

13.2724 

12.6035 

.651x10", 

.259x10", 

11.9894 

11.9954 

.131X10"° 

.557x10" 

10 

12.3605 

.112x10"^ 

12 

12.2367 

.558x10"^ 

16 

12.1228 

.194xl0"4 

TRUE  VALUE 

12.0 

12.0 

START  UP  8.0 
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Table  4.8 


(treating  all  three  parameters  as  unknown) 


N 

-N 

qi 

Linear  Spl 

q^io-3) 

ine 

— N 

q3 

J(q*) 

3 

2.25943 

.313872 

1.59287 

.181x10 

4 

4.16515 

.37462 

.83466 

.105x10 

6 

.95858 

.12019 

1.09886 

.565x10 

8 

.74126 

.13653 

1.67134 

.283x10 

12 

1.062306 

.360775 

3.17335 

.697x10 

16 

1.03356 

.504339 

4.68815 

.209x10 

20 

1.04769 

.652659 

6.08683 

.108x10 

24 

1.02567 

.745678 

7.13885 

.442x10 

28 

1.01352 

.820307 

7.98051 

.241x10 

TRUE  VALUE 

1.0 

1.2  x  lO-3 

12.0 

START  UP 

.8 

1.0  x 103 

8.0 

-3 

-3 

-4 

-4 

-5 

-5 

-5 

-6 

-6 


N 

Ciinc  Spl 
q^(*103) 

ine 

—N 

q3 

2 

2.62212 

1.69692 

5.97524 

.879x10 

3 

2.31452 

1.80666 

7.23897 

.212x10 

4 

1.00842 

1.20156 

11.91970 

.963x10 

TRUE  VALUE 

1.0 

1.2  x 103 

12.0 

START  UP 

.8 

1.0  x 103 

8.0 

The  machine  time  and  computations  required  in  use  of  our  algorithms 
are  difficult  to  give  precisely  because  this  depends  on  many  factors: 
starting  values  (initial  guesses  for  the  q^),  convergence  criteria  for 
the  optimization  method,  and  error  tolerances  in  the  ODE  solver.  Since 
our  primary  purpose  has  been  to  examine  convergence  properties,  the  con¬ 
vergence  criteria  employed  were  very  stringent  (for  example,  convergence 
to  S  decimal  digit  accuracy  in  q^)  and  thus  our  computations  were  not 


made  with  the  greatest  regard  for  efficiency.  However,  it  is  clear  that 
even  our  less-than-optimal  programs  can  be  fine-tuned  to  give  more  effici¬ 
ent  solutions  when  less  accuracy  is  needed. 

To  give  some  measure  of  the  amount  of  computation  required  and  as  a 
basis  of  comparison  between  linear  and  cubic  spline  schemes  as  used  in 
the  examples  above,  we  note  the  typical  times  required  (on  the  IBM  370 
at  Brown  University)  for  each  evaluation  of  J(q,y  )  in  the  above  examples 
ranged  from  1.5  CPU  seconds  (N  =  3)  to  11  seconds  (N  =  16)  for  the  linear 
splines  while  ranges  for  the  cubic  scheme  were  2.5  seconds  (N  =  1)  to  7.5 
seconds  (N  =  6) . 

To  compare  the  relative  efficiency  of  the  linear  versus  cubic  splines 
and  to  provide  a  further  measure  of  the  time  involved,  we  list  the  follow¬ 
ing  results  for  Example  4,  with  data  of  type  (A).  As  we  have  noted,  the 
Levenberg-Marquardt  algorithm  was  used  to  perform  the  optimization  for 
a  given  approximation  level  N.  It  is  an  iterative  procedure,  and  each 
iteration  may  require  many  evaluations  of  J  to  evaluate  the  gradient 
and  to  test  if  J  has  been  decreased  after  each  new  iterative  step.  We 
report  the  number  of  iterative  steps,  evaluations  of  J,  seconds  of  CPU 
time  per  evaluation  of  J  and  total  CPU  time  to  find  q^  from  the  start¬ 
ing  value  reported  in  Tables  4.1  -  4.3.  These  results  are  pessimistic  in 
the  sense  that  all  convergence  parameters  were  set  to  obtain  maximal  ac¬ 
curacy,  resulting  in  more  iterative  steps.  However,  it  does  serve  as 
something  of  a  basis  for  comparison.  When  treating  q^  as  unknown,  the 
following  CPU  time  was  required  (see  Table  4.1). 
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iterative 

evaluations 

CPU  seconds 

total  CPU 

approx 

N 

steps 

of  J (q) 

per  evaluation 

seconds 

Linear 

3 

5 

13 

1.8 

25.7 

4 

5 

13 

2.2 

31.4 

5 

3 

23 

2.7 

64.4 

6 

3 

9 

3.2 

29.7 

8 _ 

4 

15 

5.0 

76.5 

Cubic 

1 

5 

IS 

2.5 

37.5 

2 

3 

9 

3.0 

27.7 

3 

3 

9 

4.0 

36.0 

When  treating  q2  as  unknown,  the  following  CPU  time  was  required 
(see  Table  4.2). 


iterative 

evaluations 

CPU  seconds 

total  CPU 

approx 

N 

steps 

of  J(q) 

per  evaluation 

seconds 

Linear 

3 

7 

30 

1.5 

49.6 

4 

4 

16 

2.2 

37.0 

6 

3 

9 

3.6 

34.4 

8 

3 

9 

4.4 

44.3 

12 

3 

9 

7.2 

66.2 

Cubic 

1 

3 

9 

2.5 

22.9 

2 

3 

9 

3.0 

27.7 

3 

3 

9 

5.0 

47.3 

When  treating 

q^  as  unknown 

i,  the  following  CPU  time  was 

required 

(see  Table 

4.3) 

• 

iterative 

evaluations 

CPU 

seconds 

total  CPU 

approx 

N 

steps 

of  J(q) 

Per 

evaluation 

seconds 

Linear 

3 

5 

13 

1.7 

23.2 

4 

4 

16 

2.4 

39.8 

5 

2 

15 

2.7 

42.0 

6 

3 

9 

3.5 

31.9 

8 

3 

9 

5.0 

46.4 

10 

3 

9 

6.0 

57.5 

12 

3 

9 

7.0 

65.4 

Cii>ic 

2 

5 

9 

3.0 

27.0 

3 

3 

9 

4.0 

36.3 

4 

3 

9 

5.0 

45.3 

In  most  cases,  the  cubic  spline  approximation  was  more  efficient  than 

linear  splines.  For  example,  when  treating  as  unknown,  the  cubic 

—N 

splines  with  N  =  2  required  27.7  seconds  of  CPU  time  to  find  q^;  to 
achieve  the  same  accuracy  with  linear  splines,  it  was  necessary  to  use 
N  =  12  which  required  66.2  seconds  of  CPU  time.  In  the  case  when  all 
three  parameters  were  treated  as  unknown,  the  relatively  greater  effici¬ 
ency  of  the  cubic  splines  is  even  more  important.  The  linear  splines 
failed  to  give  a  good  approximation  even  for  N  =  16.  Some  typical  times 
to  solve  this  problem  when  all  3  parameters  were  treated  as  unknown  are 
(Table  4.4): 


approx 

N 

iterative 

steps 

evaluation  s 
of  J(q) 

time  per 
evaluation 

total 

time 

Linear 

8 

9 

31 

4.0 

126. 

16 

7 

27 

11.0 

297. 

Cubic 

8 

5 

21 

12.0 

258. 

In  summary,  both  the  cubic  spline  and  linear  spline  approximation 
schemes  converged  rapidly  when  many  spatial  observations  (i.e.  data  of 
type  (A))  were  available  and  both  gave  good  estimates  even  for  small 
values  of  N.  For  high  accuracy  estimates,  the  cubic  spline  method  was 
more  efficient.  For  problems  with  limited  data  (type  (B)  of  Example  4), 
the  difference  in  performance  of  the  methods  is  more  marked.  The  linear 
spline  estimates  for  all  three  parameters  simultaneously  appeared  to 
be  converging,  but  failed  to  give  good  estimates  of  q^  and  q^  even 
for  N  =  28;  the  cubic  approximation  on  the  other  hand  gave  a  very  accur¬ 
ate  estimate  for  N  =  4  in  far  less  time. 
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Thus  the  linear  approximations  are  satisfactory  when  high  accuracy  is 
not  required  and  when  data  is  sufficient  to  permit  linear  approximations. 
Otherwise  the  cubic  spline  approximations  provide  an  excellent  alterna¬ 
tive  which  we  believe  is  computationally  feasible. 

Conclusions 

In  this  paper  we  have  presented  parameter  estimation  methods  that 
should  prove  useful  in  modeling  large  complex  structures  as  continua. 

We  have  focused  on  systems  described  by  the  Timoshenko  theory  although 
models  for  shear  beams  (see  [1])  or  beams  with  pure  bending  (the  Euler- 
Bemoulli  theory)  can  also  be  easily  treated  with  our  techniques  (see  [4]). 
Our  methods  have  a  sound  theoretical  foundation  (i.e.  we  can  prove  con¬ 
vergence)  and  we  have  successfully  tested  them  numerically  on  a  large 
number  of  examples  (some  of  our  findings  being  reported  above) . 

We  are  currently  developing  the  methods  (both  theoretical  investiga¬ 
tions  and  computational  tests  offer  encouragement)  for  elastic  models 
with  spatially  varying  parameters  (coefficients).  We  are  optimistic  that 
these  methods  will  prove  extremely  useful  in  engineering  studies  of  large 
structures  that  are  currently  of  interest  in  aeronautics  and  astronautics. 
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